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ABSTRACT 


A computer program called HOPI was developed to predict reorientation 
flow dynamics, wherein liquid moves from one end of a closed, partially 
filled, rigid container to the other end under the influence of container 
acceleration. The program uses the Simplified Marker and Cell (SMAC) 
numerical technique and, using explicit finite-differencing, solves the 
Navier-Stokes equations for an incompressible viscous fluid. The effects 
of turbulence are also simulated in the program . HOPI can consider 
curved as well as straight walled boundaries. Both free-surface and 
confined flows can be calculated. The program was used to simulate 
five liquid reorientation cases. Three of these cases simulated actual 
NASA LeRC drop tower test conditions while two cases simulated full- 
scale Centaur tank conditions. 

It was concluded that while HOPI can be used to analytically determine the 
fluid motion in a typical settling problem, there is a current need to 
optimize HOPI. This includes both reducing the computer usage time 
and also reducing the core storage required for a given size problem. 
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SUMMARY 


This report covers the work performed under NASA/LeRC Contract NAS3-14361 by 
the Convair Aerospace division of General Dynamics during the period from February 
1971 thru April 1972. 

The objective of this contract was the development of a numerical method to predict 
reorientation or settling flow dynamics , wherein liquid moves from one end of a 
closed, partially filled, rigid container to the other end under the influence of con- 
tainer acceleration. 

Convair used the Navier -Stokes and the continuity equations for incompressible, 
viscous fluid as the basic equations governing the reorientation flow dynamics. The 
equations were programmed using explicit finite differencing for two-dimensional 
planer and three-dimensional axisymmetric problems. The Simplified Marker and 
Cell (SMAC) numerical technique was used to determine liquid-vapor interface 
positions. 

Convair usee a turbulent kinematic viscosity relationship to analytically approximate 
the effects of turbulence 

W tubb *‘ 2 Hi!t|. llrl) 

where TURB is an empirical factor which is used to correlate the data. A value of 
0.05 was found to best correlate the LeRC drop tower results. 



where DR and DZ are the grid dimensions and u and v are the velocity components in 
the radial (r) and axial (z) directions, respectively. 

The viscosity of a fluid was treated as the sum of the molecular viscosity and turbulent 
viscosity. For the reorientation flow cases run during the contract the magnitude of the 
turbulent viscosity was generally at least an order of magnitude greater than the 
molecular viscosity. 

Convair developed a computer code called HOPI using the equations formulated in the 
analytical studies. It is noted that while HOPI is restricted to two dimensional 



problems the basic SMAC technique may be three dimensional. HOPI can handle 
curved as well as straight-walled boundaries. It has the ability to calculate both 
free -surface and confined flows. It can be used in either cylindrical or plane geometry. 
The size of the computing mesh is easily changed from problem to problem. The grid 
dimension in each direction must be constant through the grid mesh, although the grid 
dimensions in the radial and axial direction may differ. Gravitational effects may be 
included in any orientation however it is noted that HOPI can only compute axisymmetric 
flow when a cylindrical geometry is used. For any given time interval the gravitational 
acceleration must be constant. HOPI has a surface pressure interpolation scheme that 
prevents unrealistic breakup of the surface. 

Convair used the computer program HOPI to simulate five liquid reorientation cases. 
Three of the cases simulated actual NASA LeRC drop tower test conditions and were 
used to provide confidence in the results generated by HOPI. The additional two cases 
simulated full-scale Centaur tank conditions. The cylindrical tanks had height to 
diameter ratio of two, hemispherical forward and spherical segment bottom dome, 
with a radius of 7 cm for drop tower cases and 150 cm for full scale cases. In all 
cases the initial interface shape corresponded to a Bond number of 10. The 
reorientation Bond numbers ranged from 100 to 450. 

The experience gained in running these cases indicates that HOPI can be used to 
analytically determine the fluid motion in a cryogen storage tank under a continuous 
settling load. However, there is a current need to optimize HOPI. This includes 
reducing the computer usage time and also reducing the core usage required for a 
given size problem. To reduce the computer usage time it is recommended that 
HOPI be modified to handle variable grid mesh and to change the subscripting from 
double to single throughout the code. To reduce the core storage it is recommended 
that overlaying be used. 
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1.0 INTRODUCTION 


Auxiliary thrustors are used to create reduced-gravity environments to control and 
locate liquid surfaces in space vehicle systems. The continual use of even small 
thrustors for long-term missions can result in undesirable weight penalties and, 
thus, their operation may have to be intermittent. During the off -times ,the liquid 
surface may be destabilized and relocated due to disturbances acting on the vehicle. 
Various passive retention devices such as capillary baffles might be used to control 
the location of the liquid, but it appears that, at least for large cryogenic propellant 
systems, reduced-gravity auxiliary thrustors will remain a primary method of 
control. The auxiliary thrustors will be relied upon to settle or reorient the pro- 
pellant back to its desired location prior to restart or venting operations. The thrust 
and time required to settle efficiently and reliably are vital design factors. The 
initial destabilization and ensuing flow dynamics have been studied analytically, but 
reorientation time estimates are still derived almost entirely from empirical results 
based on scale-model experiments. References 1 through 7 discuss aspects of this 
reorientation problem. 

During this contract, a computer program called HOPI was developed which can be 
used to predict reorientation flow dynamics, wherein liquid moves from one end of 
a closed partially filled, rigid container to the other end under the influence of 
container acceleration. This computer program numerically solves the Navier- 
Stokes equations for viscous incompressible fluids using the Simplified Marker-and- 
Cell (SMAC) technique. The basic development of this technique is given in Reference 
8. One limitation of the program described in Reference 8 is that it cannot consider 
curved boundaries . Basic techniques to handle curved boundaries were developed by 
Viecelli in Reference 9. 

The objective of this study is to extend the scope of past analytic studies of 
reorientation to the liquid-liquid impingement phase and to calculate during this phase 
interface profiles and liquid volume collection rates. The region of particular interest 
is reorientation Bond numbers of 100 to 500. Representative cases were calculated 
to compare directly with drop tower experiments and to corroborate empirical results 
for full-scale vehicle propellant tank configurations. 
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2.0 DISCUSSION 


In the Simplified MAC (SMAC) technique, the procedure for a calculational cycle is as 
follows : 

1 . A tentative field of advanced-time velocities is calculated by using an arbitrary 
pressure field within the field, but with a pressure boundary condition at the free 
surface satisfying the normal stress condition. Correct velocity boundary conditions 
assure that this tentative velocity field contains the correct vorticity at every interior 
point in the field. The tentative velocities do not satisfy continuity (i.e. , V • V = 0). 

2. The tentative velocities are modified to their final values so as to preserve the 
vorticity at every point. A potential function is employed, determined by the 
requirement that it convert the velocity field to one which satisfies the incompress - 
bility condition everywhere. 

The computer program developed during this contract is called HOPI, Reference 10. 
It handles two-dimensional or axisymmetric three-dimensional problems involving 
incompressible Newtonian fluids; In the following sections the basic equations and 
techniques used in HOPI are presented. 

2.1 EQUATIONS USED IN HOPI 

The basic differential equations are 

du + J_ ar a u 2 + duv 
at tP ar 9Z 




d0 

d Z 







The velocity components, u and v, are respectively in the r and z directions, and the 
pressure, 0 , is normalized to unity density. The two components of gravitational 
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acceleration are g r and g z as indicated below. Plane (Cartesian) coordinates have 
a = 0, and in cylindrical coordinates, a - 1. 0. 


z 



In differential form, the equation for transport of vorticity, u) , is independent 
of the pressure, so that any field of pressure inserted into the Navier-Stokes 
equations will assure that the resulting velocity field carries the correct 
vorticity. An arbitrary pressure field will not, however, assure the vanishing 
of D, but if the velocity field is altered by the addition of the gradient of an 
appropriate potential function, the resulting field will carry the same vorticity, 
have vanishing D, and accordingly will be uniquely determined, hence correct. 

This is the essence of the procedure for the finite difference solutions. As a 
starting point, the equations are written in the following form. 


u i+l/2,j ~ u i+l/2, j = 
6t 


r a„n , „n _or „n ,.n , 

r i u i+ l/2, j u i-l/2,i ~ r i+l u i+3/2. j u i+l/2, j 


r Sl/2 6r 


n n n n 

U i-fl/2,j-l/2 V j+l/2, i-l/2 ~ Vl/2, j+l/2 Vl/2, j-H/2 


6z 


- ib 

+ Id HlJ + g r 

6r 


+ v 


1 , n 


n 


n 


6z ? ^ U i+l/2, j+1 + U i+l/2,j-l 2U i+l/2,j ^ 


1 n n n 

6r6z V i+1, j+l/2 V i+l,j-l/2 V i, j+l/2 




( 1 ) 



and 


~n n 

V j,j+l/2 V i, j+1/2 

6 t 


an n an n 

*1-1/2 “i-l/2 J+1/2 V i-l/2 J+1/2 ~ r i+l/2 “ i+ 1/2, j+1/2 V j+l/2, j+1/2 

r . or 
1 


n n 

v. . _ v 


n n 

- v . . _ v . 


. i J+1/2 i, ,1-1/2 U+3/2 i, .1+1/2 

6 z 


+ . \ l ~ 

6 z 


& z 


n 
u . 


n 


V— f r a , “i+1/2, j+1 
a 6 r L i+1/2 6 z 


- u 


1+1/2J 


r. 6 r 
i 


) 


n n 

V i+1, j+1/2 ~ V j J+1/2 
6 r 

n n 

a i-l/2,]+l 1 - 1 / 2 , j 

i-l/2 1 6 z 


n n 

V i, j+1/2 ~ V i-l J+1/2 
6r 


( 2 ) 


D 


a n+1 a n+1 

n+1 _ r i+l/2 U i+1/2J r j-l/2 U i— 1/2, j 
ij 

r, 6 r 
i 


n+1 n+1 

+ V i J+1/2 ~ V i,j-l/2 
6 z 


= 0 


<3) 
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The subscripts refer to position in the finite -difference mesh (see Fig. 1), and the 
superscript n counts time cycles. The true, pressure, 0 , has been replaced by 
the arbitrary field, 0 , and accordingly the new-time velocities are marked with 
tildes . 


j+ 1/2 

j 

( « 

j-1/2 


CELL (ij) 


V ij+l/2 



Figure 1. The Location of the Cell Variables in a 
SMAC Cell 

The finite difference equations have been written with explicit (retarded-time) fluxes . 
Cell-centered momentum convection terms have been written in the ZIP form, which, 
although continuing to assure internal momentum conservation as in MAC, allows 
SMAC to conserve momentum in the immediate vicinity of a rigid wall. This type 
of differencing introduces the definition 



U i-l/2,j U i+l/2,j 


(4) 


An advantage of its usage is the removal of a destabilizing truncation error term 
that occurs in the original form of MAC, (Reference 11). 


A finite-difference approximation to the vorticity is 


n 

°i+ 1/2, j+ 1/2 s 


n n 

U i+l/2, i+1 U i+l/2,j 
6 z 


n n 

V jH, j+l/2 " V i, i+l/2 


dr 


(5) 
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with centering at cell corners. Equations (1) and (2) can be combined to obtain 
a transport expression for <u. + jy 2 j+i /2 w ^ ic ^> like the differential equation, is 

independent of the 4> field. Accordingly, the explicit calculations of the tilde 
velocities assures that the vorticity at ^every internal mesh corner point is correct, 
independent of the choice of ip ■ This is not true, however, for corner points that 
lie on rigid walls, which are not correct until the tilde velocities have been corrected 
to assure the vanishing of D. For purely explicit calculations, which are acceptable 
for Reynolds numbers greater than about unity, SMAC vorticity diffusion from the 
wall is nevertheless correct, because the tilde velocities are based entirely upon 
the final velocities from the previous cycle, which do agree with the proper wall 
vorticity. 


A cell is flagged as a surface (SUR) cell when it contains fluid, marker particles 
and it has at least one adjacent neighboring cell which is flagged empty. Marker 
particles do not perform any function in HOPI calculations other than to indicate 
the position of any free surface that may be present. On free surfaces the 
tangential stress condition is 


3u 5 v 
3 z 3r 



so that u i+1 /2 j +1 is determined by the equation 


u 


= u 


dz 


i+l/2,j+l “i+l/2, j <5 r ( Vl, j+l/2 " V i, j+l/2 


- v. 


,)■ 


( 6 ) 


This assures that the tangential viscous momentum flux vanishes when calculated 


by Equation (1) for u. 


i+1/2, j 


Iu addition, the normal stress condition is 




(applied) + 


2 v 

dz ^ V i, j+l/2 


v 

i, 


j-l/2 ) 


( 7 ) 


The applied part of the pressure is specified according to the requirement of the problem 
while the viscous part assures that there is otherwise no net flux of normal momentum 
through the surface. It is important that the normal stress condition be placed on the 
free surface rather than at the center of the surface cell. 
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The viscosity coefficient u used in Equations 1 and 2 is the sum of the kinematic 
molecular viscosity and the turbulent viscosity. 

u ^molecular + u turb 

The molecular viscosity is an input quantity and is a fluid property. The turbulent 
viscosity coefficient is calculated internally in the program as indicated below. 



and TURB is an input quantity. 

This expression for turbulent viscosity is of the form predicted by both Prandtl's 
mixing-length theory and Taylor's vorticity transport theory. While other expressions 
for predicting the turbulent viscosity do exist, the above was selected due to its wide 
acceptance and simplicity. 

Basically a turbulent viscosity is calculated in a cell containing fluid when at least 
two of its adjacent neighboring cells also contain fluid. The criteria of requiring 
fluid in adjacent fluid cells is needed so that dv/dr and du/dz can be calculated. 

For the reorientation flow cases run during this contract (Section 3) the magnitude 
of the local turbulent viscosity coefficient was at least an older of magnitude greater 
than the molecular viscosity coefficient for most of the duration of the problem. This 
indicates that the viscosity coefficient used in Equations 1 and 2 is mainly a result 
of the turbulent viscosity coefficient. 

2.2 HOPI COMPUTER PROGRAM 

This section describes in detail the SMAC calculational cycle in the framework of 
HOPI. HOPI embodies a number of features that make it a useful tool. Among these 
are 

1. It is written in FORTRAN IV for the CDC-6400 computer. 

2. It has the ability to calculate both free-surface and confined flows. 

3. It has the ability to calculate either in cylindrical or plane geometry. 
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4. The size of the computing mesh is easily changed from problem to problem. 

5. It has a simple, straightforward setup, allowing different initial conditions and 
particle resolution in different regions of the mesh. 

6. Various boundary conditions are available, along with an obstacle. 

7. Gravitational effects may be included in any orientation. However, it is noted that 
HOPI can only compute axisymmetric flow when a cylindrical geometry is used. 

8. Both curved and straight wall boundaries can be used. 

The underlying HOPI scheme is now discussed in detail. 


2. 2. 1 COMPUTING MESH. The HOPI computing mesh will handle either cylindrical or 
plane geometry calculations. The equations are all presented in cylindrical form, however, 
it is shown in Section 2.2. 6 that it is simple to transform them to plane form. Hie 
computing mesh as used in cylindrical coordinates (Figure 2) is an infinitely thin radial 

radial slice from a cylinder. The computing 
cells become "toroids" of revolution about 
this cylinder, and are of uniform size in 
I CELLS HOPI. 







i 

__L 




6z 




T 




-* 

6r 

•*— 










1 


J CELLS 


The radial coordinate is denoted by r, and 
the radial cell size is denoted by 6r, 
whereas the axial coordinate is denoted by 
z, with 6z the cell size in the axial 
direction. The origin is located at the 
lower left corner of the mesh of cells, 
the centers of which are labeled with the 
indices i and j. These indices increase in 
the r and z directions, respectively. The 
mesh of cells isT (IBAR) cells wide and J 
(JBAR) cells high. 


-► r 


Figure 2. Computing Mesh in Cylindrical 
Coordinates 


2.2.2 CELL VARIABLES . The primary 
cell variables are the two components of 
velocity, u and v, the pseudo-pressure K 
and the velocity divergence, D. The 
centering of these variables about a SMAC 
cell is shown in Figure 1, where it is seen 

that the velocities in the radial direction are centered at the left and right edges of the 
cell, and axial velocities at the top and bottom edges. D, and tjr are defined at the cell 
center. In the equations, velocities appear that are not positioned at the cell edges. In 
these instances, simple averages are calculated. For example 
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(uv) . 

1+1/2, J— 1/2 


BASIC TYPES OF CELLS AND THE INDEXING SCHEME. 
cells which may exist within the mesh are described in Figure 3. 
is demonstrated in Figures 4 and 6. 


The basic types of 
Use of these flags 


ARB An OB cell which has a fluid particle within e • DR of the arbitrary boundary, where e is the boundary 
sensing parameter and is usually set equal to .25. 

BOT A COR cell containing a segment of an arbitrary boundary which has its midpoint and an$e being stored in 
the OB cell just below it. 

COR, A cell which has a line segment of the arbitrary boundary passing through it, however fluid area to total cell 
area fraction is less than .25. The fluid area is to the left of the line segment. Also, any cell just outside 
an OB cell is a COR cell. 

EMP The cell is empty (contains no fluid particles). 

EOC A cell which is either EMP or COR. 

EXT Any cell outside a COR cell. 

FUL . A cell which contains fluid and has empty neighbor. 

LEF A COR cell containing a segment of 
an arbitrary boundary which has its 
midpoint and angle being stored in 
the OB cell just to the left of it. 


OB A cell which has a line segment of 
the arbitrary boundary associated 
with it (see COR). 


OK An ARB cell that does not contain 
the intersection of a free surface 
and the boundary. 


RIG A COR cell containing a segment of 
an arbitrary boundary which has its 
midpoint and angle being stored in 
the OB cell just to the right of it; 

SNC A cell which is also flagged as FUL 
or SUR and not COR. 


SUR A cell which contains fluid and has at 
least one empty neighbor. 

TOP A COR cell containing a segment of 
an arbitrary boundary which has its 
midpoint and angle being stored in 
the OB cell just above it. 


EXT 

EXT 

EXT 

EXT 

EXT 

EXT 

EXT 

B0T 

B0T 

'-'"ut, 

B0T 

t-U-Lu*.. 

C0R 
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0K 

are/ 

0k 


W 

Arb, 

0K 

FUL 

FUL 

FUL 

FUL 

FUL 

FUL 

FUL 

SUR 

SUR 

FUL 

FUL 

FUL 

FUL 

FUL 

EMP 

EMP 

SUR 

SUR~^ 

SUR 

FUL 

FUL 

EMP 

EMP 

EMP 

EMP 

EMP 

SUR^ 

SUR/ 

EMP 

EMP 

EMP 

EMP 

EMP 

EMP 

EMP 

0B 

0B 

EMP 

EMP 

0B 

0B 

0B, 

-rrrnrrr 

C0R 

T0P 

~^yg~ rrn 

rrrrrrrr 

0B 

"tHp* 7 

T0P 

C0R 

EXT 

EXT 

C0R 

C0R 

EXT 

EXT 

EXT 


-Arbitrary 

Boundary 

Region 
Containing 
f Liquid 

-Liquid- 

Vapor 

Interface 


} 


Region 

Containing 

Vapor 
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NOTE: 

1. All ARB cells are also flagged as 0B. 

2. All T0P, B0T.RIG and LEF cells are also flagged as C0R cells. 

3. All the above SUR and FUL cells are also flagged SNC. 


Figure 3. Grid Network With Typical Cell Flags 




In Eulerian programs of the MAC-SMAC type, experience has shown that it is 
helpful to surround the above-described mesh with a belt of cells lying outside 
the boundaries. Cells in this belt are called boundary (BND) cells, and they 
simplify the task of handling boundary conditions, particularly those relating 
to velocities. 

In the FORTRAN code HOPI involving the SMAC equations, we reference "'I'i, j" 
simply as "PSI (I, J). " But "Ui+i/ 2 , j" cannot be referenced by a."half-integer" 
index in FORTRAN, so the convention has evolved that "U(I, J)" refers to this 
velocity, "V(I,J)" in the code actually refers to j+i/ 2 * 

Given this referencing scheme, the reason becomes apparent for creating a column 
of BND cells on the left for storing the u's (normal velocities) for the left boundary , 
and a row of BND cells along the bottom for storing the v's (normal velocities) 
pertaining to the bottom boundary. In addition, tangential velocities are logically 
arranged in BND cells on all four sides; their values are based on neighboring u's and v's 
lying inside, modified through the existing wall boundary condition. With appropriate 
"outside" velocities thus distributed over the BND cells, the computer code can solve 
the SMAC equations over the interior cells, picking up velocities on all sides as 
needed, Without having to test for adjacency of a mesh boundary. Variables in BND 
cells must, however, be updated as the corresponding variables inside change in value. 

It is clear that it is possible for a cell to contain a number of flags. For example, 
a cell could be OB, ARB, FUL, OK, and SNC all at once. In HOPI an NBIT function 
is used to determine if a flag is set for a given cell. To speed up the computation 
time certain cell flags such as EOC and SNC were developed, which represent two or more 
flags. To further speed up the computation time "G" flags were developed as indicated 
below: 

G = 2 implies a COR cell 

G = 3 implies an OK cell ' 

G = 4 implies EMP and not COR or EXT cell 
G = 5 implies a BND cell 

An optimization study indicated that each computation of the form 
IF (NBIT (10, F(N)) . EQ. 1) GO TO 100 
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takes . 386/10, 000 seconds while a computation of the form 

IF (G(N) . EQ . 1) GO TO 100 
takes . 141/10, 000 seconds. 

This shows that it is 2. 75 times faster if a G cell flag is used instead of the NBIT 
function. A given cell can contain only one G flag. 

2.2.4 MARKER PARTICLES . In addition to the mesh of Eulerian cells, SMAC 
employs a set of massless marker particles, which are helpful for allowing a visual 
representation of the fluid, but whose essential purpose is to define the position of 
the free surface so that the configuration of SUR cells can be sensed. Beyond this, 
the marker particles do not enter into the calculation, but are merely embedded in 
the fluid and are carried along by it. Each cycle the marker particles are moved 
with a weighted average of the four nearest u's and of the four nearest v's. 

2.2.5 BOUNDARY CONDITIONS . In HOPI both straight line and curved wall 
boundaries can be used. A curved wall is approximated by a series of straight line 
segments within the grid mesh, where each cell that has part of the curved boundary 
passing through it contains a straight line segment. Each line segment is formed by 
joining the two points formed where the curved boundary crossed the rectilinear 
Eulerian boundary of the cell. 

Straight Wall Adjacent to a BNP Cell 

With reference to Figure 4, the indices i, j refer to the cell inside the system, and 
i-1, j refer to the BND cell lying just outside. 


OUTSIDE FLUID SIDE 
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1. FREESLIP: A freeslip boundary represents an axial centerline or a plane 
of symmetry or a non-adhering surface that exerts no drag upon the fluid. 

The normal velocity component vanishes at the wall, and there is no gradient 
in either tangential velocity or in the potential function ip 


" U i-l/2,j =0 

Vl,j+l/2 = V i,j+l/2, 

I V i-1, j-1/2 =V i,j-l/2 




2. NOSLIP; A noslip boundary represents a viscous boundary that exerts a drag 
upon the fluid. This is accomplished by forcing the tangential velocity to go to 
zero at the wall 


Vl/2,j ' 0 


V i-l,j+l/2 V i,j+l/2 , 


Vl, j-1/2 _V i, j-1/2 , 

i =4 / 

i-l.J M • 



The no slip option is not available for curved wall boundaries. Also, it is noted that 
often more accurate results might be obtained by imposing the free slip condition on a 
boundary even when in reality a no slip condition exists. This is true when the grid 
mesh is too coarse to accurately resolve the detailed motion of the boundary layer 
resulting from a no slip condition. This was true during the cases run during this 
contract. Therefore, for the cases reported in this report the free slip boundary 
condition was always imposed on every boundary. 

Curved Wall Boundaries 

The basic technique for handling curved-wall boundaries was developed by Viecelli in 
Reference 9. Only the case of free-slip boundary conditions exist for curved- wall 
boundaries. The motion of an interface between a liquid and a curved wall 
is equivalent to that of a free surface with an applied pressure distribution. Given 
any interfacial shape or motion one can produce that same shape or motion with some 
unique pressure distribution applied to a free surface. This is the basis for arbi- 
trary boundary calculations . 
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A curved wall is specified by a locus of points. These points might be the edge 
intersection points of a mesh covering the wall or body; however, in general the 
points will not lie on any of the Eulerian mesh lines used in the finite difference 
solution of the hydrodynamic equations. The first step is to connect successive 
points with straight lines and find all of the intersections of these segments with the 
underlying Eulerian mesh. One can then represent the boundary section lying inside 
an Eulerian cell by a single straight line connecting the points of intersection of the 
boundary with the sides of the cell. It is possible to conceive of situations where the 
boundary winds in and out of a single cell creating some ambiguity. However, this 
occurs only when there are an insufficient number of zones to accurately calculate 
details of the flow. Therefore it is assumed that the Eulerian zoning is always fine 
enough so that the boundary has only two intersections wife each cell. Having broken 
the boundary into a set of straight line segments, each associated with a unique 
Eulerian cell, one can specify the position of each segment by a unit vector normal to 
the boundary segment positioned at the midpoint of fee segment. HO PI uses the 
convention that the normal points towards the liquid and to the left as one advances 
from the ith to i + 1th boundary point. 


The second step is to define and flag Eulerian boundary cells as those along the inside 
edges of the contour approximated by the boundary segments . The liquid area of the 
boundary cells are then calculated again using the convention that the liquid is to the 
left as one advances from the ife to i + 1th boundary point. Then, if the liquid 
fraction of the total cell area is greater than 1/4, the boundary cell flag, OB, is 
turned on. If the liquid area fraction is too small the cell is flagged COR. It is then 
determined to which of the four neighboring cells the boundary segment normal points 
nearest and then set the boundary flag, OB, for that cell. When that cell also contains 
a boundary segment the two segments can be replaced with one by removing the 



Figure 5. Replacing Two Boundary 
Segments With One Segment 
Spanning Two Cells 


boundary intersection point between the two adjoin- 
ing cells. This new segment, spanning two cells, 
is also defined by velocity and position vectors at 
its midpoint. ThiB case is illustrated in Figure 5. 
Note that as a result of this selection process the 
midpoint of a boundary segment may not physically 
lie in the Eulerian boundary cell with which it is 
associated. Thus, it is sometimes necessary to 
assign a pointer to a boundary cell indicating which 
neighboring cell contains the midpoint of the 
associated boundary segment. Once the boundary 
cell flags have been set it is necessary to go 
through them and turn off the boundary (OB) flag 
in any cell in the corner of a right angle cell 
pattern, since such cells are bounded on two 
adjacent sides by either a pair of interior or 
exterior cells and on the other two by boundary 
cells. Figure 6 shows a closed boundary curve, the 
approximate Eulerian boundary, and the cell pattern. 
It illustrates the various cases mentioned above. 
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This process of defining a set of Eulerian boundary cells is analogous to the defining of 
a set of free surface cells in the MAG method but with the additional constraint that the 
resulting cell pattern avoids overdetermining the boundary condition on the pressure. 
Cells just outside the Eulerian boundary are flagged as COR cells while the remaining 
cells outside the Eulerian boundary are flagged as external (EXT) cells. This is so that 
one can tell when velocity components lie on the exterior sides of boundary cells. 
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NOTE 

A curved wall is approxi- 
mated by a series of 
straight line segments 
within the grid mesh , where 
each cell that has part of 
the curved boundary passing 
through it contains a 
straight line segment. 

Each line segment is 
formed by joining the two 
points formed where the 
curved boundary crossed 
the rectilinear Eulerian 
boundary of the cell. 










Figure 6. Permanent Cell Flags for Arbitrary Boundary 
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Once the boundary cells (flagged - OK) have been determined, the following 
relaxation equation is used to compute the pressure in these cells: 



( 8 ) 


In this equatioji n is the normal defining the boundary segment associated with cell 
(ij) and (\£ +1 ) k is the liquid velocity at the midpoint of the sgement computed with the 
MAC area weighted interpolation formula as shown in Figure 7. 

Clearly (vj| ). j is one of the iterates and must be recomputed each time the 
pressures and velocities are adjusted. The relaxation parameter and minimum 
mesh dimension are RELAX and X respectively. The formula shows that instead 
of adjusting the pressure proportional to the divergence or net flux out of a cell 
one adjusts it proportional to the flux across the boundary measured relative 
to coordinates fixed in the boundary. If liquid is flowing across the boundary the 
pressure will be increased until the outflow stops. Conversely, if liquid is tending 
to separate from the boundary the pressure will decrease until the liquid flows 
tangent to the boundary. 


Velocities at the exterior sides of the Eulerian boundary cells and other exterior 
points are necessary in the area weighting formula, and must be recomputed during 
each iteration sweep. These velocity components are determined in the same way 
as those on the open sides of free surface cells in the MAC method. In Figure 8 
representative sections of boundary are shown in more detail. The original Lagrangian 
boundaries are indicated by dotted lines, and the resulting boundary segments and normals 
are shown. Heavy dark lines outline the outer edges of the boundary cells. The 
values of the velocity components at the edges and outside the Eulerian boundary, 
where necessary, are given in terms of their interior neighbors for the 2-dimensional 
plane case. In 3-dimensions with axial symmetry radius factors would be necessary 
to preserve continuity. 

In addition to calculating new cell pressures during each cycle of iteration one must 
also recalculate the velocity components. During each iteration the sum of the old 
velocity component at time nAt, the advection and the viscous terms are stored in 
u i+ l/ 2 , j and j+ 1/2 which need to be computed only once. Changes in the new cell 
velocity iterates then depend only on changes in the gradient of the pressure iterates. . 

In HOPI the pressures of arbitrary boundary cells are located at the boundary. 

This requires that one interpolate to find the cell centered pressure to use in the 
momentum equation. Figure 9 illustrates the method of linear interpolation. A 
real advantage of the simultaneous iteration method and Equation 8 is that it gives 
a simple automatic way of including complicated boundary conditions. 
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Figure 7 . Interpolating Scheme for Obtaining the Liquid 
Velocity at the Midpoint of A Boundary Segment 
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Figure 8. Determination of Mesh Velocities on the Exterior Sides of Eulerian 
Boundary Cells and Assignment of Boundary Segment Normals 



Z 


-► r 


♦a = (a * 4 + DZ ^)/(a + DZ) 

* B = *2 + b( VV /(DZ - b) 

♦ C = (c tg + DZ t 3 )/(C + DZ) 


Figure 9. Pressure Interpolating Scheme 



The marker particles, which are typically input at a density of atjeast four per cell, 
specify the fluid configuration with an uncertainty much less than the Eulerian mesh 
width. Because of this some finer criteria other than just the knowledge that a boundary 
cell contains particles is necessary. We require in addition that 

(Xp - -V • 

where Xp is the particle position, X n is the position of the midpoint of the boundary 
normal, 'n'is the boundary normal, and e is some fraction of the cell width \ typically 
1/4. Thus, we do not begin computing a pressure in boundary cells until the particles 
come within e \ of the boundary segment. 

When free surfaces are present we also need to know how to treat cells containing 
the intersection of curved wall boundaries and free surfaces. The pressure at the 
intersection point should be equal to the ambient pressure, but because the pressure 
is defined only on the Eulerian net, it is sometimes not possible to zero the flux at 
the boundary consistent with vanishing divergence without introducing a pressure. 

This happens when the angle between the free surface and the boundary is small and ; 
the liquid is colliding with a wall, producing a liquid layer on a scale too fine to be 
resolved by the Eulerian mesh. We define an intersection cell to be one that contains 
liquid and has one or more empty interior or pressure surface neighbors, and one or 
more exterior neighbors i When this definition is satisfied, the pressure is set equal 
to ambient pressure and the velocities are adjusted directly. In most circumstances 
the liquid in the cell will be part of a much larger mass. When there are one or two 
liquid neighbors, the velocity components at the sides in contact with the liquid are 
preserved, and those at the open and boundary sides adjusted to make the velocity 
tangent at the boundary consistent with vanishing divergence. In the case of one liquid 
neighbor, the velocities at the opposite cell sides are assumed equal, and the compo- 
nent with both sides open or boundary is adjusted. In either case the flux at the boun- 
dary is a linear function of a single variable, and the zero is easily found. If the 
velocity at the boundary is initially directed away from the boundary, nothing need be 
done. The remaining possibility is that there are no liquid neighbors, as happens 
when a small isolated element strikes the boundary. In this case we set the component 
of the particle velocity normal to the boundary equal to zero, and preserve the tangen- 
tial component. If a gravitational force is present we accelerate the particle velocities 
by the component of the gravitational vector tangent to the boundary. This is a free 
slip condition. 



2.2.6 CALCULATIONS TO PLANE COORDINATES . The equations in HOPI are all 
written in cylindrical coordinates, but the program is arranged so that the calculations 
can be performed in plane coordinates with a minimal loss in computing efficiency. 

The geometry type is specified by the user through the input quantity PC , which equals 
0. 0 for calculations in cylindrical coordinates, or 1. 0 for calculations in plane coordi- 
nates. For plane coordinates, fix is input in place of fir and fiy in place of 6 z. 
From this point on, HOPI performs the transformation automatically in the following 
manner. 

The equations contain coefficients involving the use of rj, the radius to the center of 
cell i, and ri±i/2 ra dius of right and left edges of cell i. To avoid recalculating 
these radii each time they are encountered, HOPI generates a number of tables which 
are entered with the index i for the quantities r i> 1 / r i+i/2» l/ r i* ( 4r i " 6r)/ 

(4r- + Sr), and (4r. +6r)/(4r i - Sr). If plane geometry is specified by PC, however, 
every entry in all six tables is generated as unity, causing all radial effects to 
disappear. Further; PC itself directly appears in some of the equations. For 
example, an expression applicable only to cylindrical coordinates is multiplied by 
(1. 0~PC), forcing its cancellation in plane coordinates. 

2.2.7 PROBLEM SETUP . Data punched on input cards provide the setup routine 
with the information required to generate the flow field at initial_time. Basic require- 
ments are the number of interior cells in both directions (I and J), the geometry 
(cylindrical or plane), the size of each cell (fir by 6 z or fix by fiy), the initial layout 
of marker particles specifying the number of particles for each cell, and the initial 
u's and v's for those cells containing particles. Other necessary information is 
included that specifies the boundary types on the four sides of the mesh, the location 
of an obstacle if one is present, the fluid viscosity, the amount of gravitational accel- 
eration, the time step (fit), and the intervals at which plots or prints are to be made. 

In the HOPI set up, the mesh is initially flagged as all EMP cells, surrounded by 
BND cells. If an obstacle is present, its cells are permanently flagged OB. Cells 
outside OB cells are permanently flagged either COR or EXT. Subsequently, any 
cells containing particles are flagged FUL. At this point, a set up containing a free 
surface has no SUR cells, but this is remedied immediately in the first ealculational 
cycle when the cell flags are adjusted. 

Finally, the setup routine must calculate scaling parameters for the microfilm plots 
that will be made of the marker particles, and calculate various coefficients and 
tabular data that will be used repeatedly in the ealculational cycle. 

At this point, the setup ends and control passes to the first ealculational cycle. Omis- 
sion of any variables, such as SUR cell flags, tangential velocities in BND cells and on 
the free surface, is rectified in the first cycle, because their calculation is a standard 
part of every cycle. 
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3. 0 " RESULTS AND DISCUSSION OF RESULTS 

During this contract five liquid reorientation cases were run using the computer 
program discussed in Section 2. Three of the cases simulated actual NASA LeRC 
drop tower test conditions and were used to give confidence in the results generated 
by HO PI. The additional two cases simulated full-scale Centaur tank conditions. 

3. 1 CORRELATION OF ANALYTICAL RESULTS WITH TEST DATA AND THEORY 

Three analytical runs were made which corresponded to actual LeRC drop tower test. 
These runs were used to help establish the accuracy and limitations of the existing 
code to solve typical liquid reorientation problems. All three problems had the 
following characteristics: 

• Height to diameter ratio of 2 

• Cylindrical tank with a radius of 7 cm 

• Fluid (FC-78) had a kinematic viscosity v = 4.74X10 -3 cm 2 /sec (Ref. 12) 

0 Initial interface shape corresponding to an initial bond number of 10 

• Hemispherical forward dome 

• Initially the fluid was at rest (all velocities = 0) 

• Grid mesh of DR = DZ = 1 as indicated in Figure 10 

• All boundaries were input with the free slip condition (Section 2.2.5) 

In all plots shown in this section just half of a tank is shown due to symmetry. The 
left side of each plot corresponds to the center of the tank. 

CASE 1 . Case 1 had the following additional characteristics: 

20% liquid (by volume) 

Spherical segment bottom dome 
Settling acceleration (GZ) of 31.4 cm/sec 2 
Turbulent viscosity coefficient = 0.0 

Initial Marker particle density =4x4 per cell (580 particles) 

Figures 11A through 11C give computer generated SC 4020 plots of the Marker particles. 
Figures 12A through 12C give computer generated SC 4020 plots of the velocity vectors 
at selected times. Initially the fluid was at rest (all velocities = 0). It is noted that 
in Figure 12A certain regions which contain fluid do not have a velocity vector indi- 
cated. This results since the velocities must be larger than a certain minimum value 
to be printed as a SC4020 velocity vector. 

The results agree closely to those of actual test data supplied by NASA LeRC. The 
liquid flows down along the side of the cylinder with a velocity V = .87 x ,g z x t 




23 













(Reference 2) where t is the time duration that the settling g-load is applied. The 
liquid in the center of the tank moves up toward the top bulkhead and the fluid adjacent 
to the top bulkhead moves parallel to the bulkhead as indicated by the velocity vectors 
in Figure 12 C. 

CASE 2 . Case 2 had the following additional characteristics: 

70% liquid 

One cell truncated spherical segment bottom dome 
Settling acceleration of 31.4 cm/sec 2 
Turbulent viscosity coefficient = . 05 

Initial Marker particle density =6x6 per cell (4285 particles) 

This problem was run for approximately 2.4 seconds of settling time. Some problems 
which did not occur in the first run appeared. The first problem was at the corner 
where the fluid coming down the wall impinged on the bottom dome. It was found 
necessary to truncate the dome by at least one cell since, for the present grid size, 
the program could not resolve the fluid motion. This is illustrated schematically 
in Figure 13. Figures 13A and 13B are the corner without and with one cell trun- 
cation respectively. 

The vector arrows indicate the motion of the fluid. In Figure 13A it is seen in cells 
8, 4 and 8, 3 that the fluid moves in two different directions. Since for a given cell 
there is only one set of velocity for all the particles in the cell, the two velocities 
cannot be resolved. 

This problem can be corrected by using a smaller grid mesh. However, since the 
present coding does not allow for variable grid sizes, the finer grid size needed 
would greatly increase the running time of the problem. 

A second problem also developed when the test case was run past the point of initial 
impact with the bottom dome. The problem was in the surface definition as defined 
by a subroutine called SURP, SURP has two main functions: one is to evaluate 
surface tension forces, the second is to extrapolate the pressures as calculated at the 
surface to find the appropriate pressure for the center of the surface cell. Figures 
14A and 14B illustrate how the surface becomes defined along the arbitrary boundary. 
This results since the particles that initially hit the bottom stagnate while the 
particles that follow move across the boundary. Because of this phenomena it is an 
inherent limitation of the program that subroutine SURP cannot be used once a 
condition like this exists. 

The first approach to solving this second problem was to delete all surface marker 
particles. However, as shown on Figure 15A this resulted in an instability in the 
surface. The second approach consists of deleting surface marker particles as they 
entered cells containing segments of an arbitrary boundary. As indicated in Figure 
15B the instability is gone. It is therefore recommended that this second approach be 
used. 
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Figure 13. Magnified Section of Grid Network Showing Cells Containing 
Arbitrary Boundary Segments 
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Figure 14. 


Figure 15. 




Surface Definition After Impact of Fluid on Arbitrary Boundary 




Effects of Surface Marker Particles in Ullage Section for Test 
Model With 30% Ullage at 1.05 Seconds 
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An examination of the geyser that has formed on Figure 15B illustrates another 
problem: a condition of rarified particles. This problem results when the particles 
are accelerating along the bottom boundary. This acceleration results when fluid is 
moving in toward the center of a cylinder under its momentum. As the radius of the 
cell containing fluid decreases, the radial component of velocity in that cell increases 
inversely with the radius. 

The condition of rarified particles becomes a problem when a cell that should contain 
particles does not. This is the case in the geyser seen in Figure 15B. A subroutine 
(called ADPART) was written to remedy this problem. Figure 16 shows a typical 
geyser with the additional particles as added by subroutine ADPART. 

Also while running Case 2 the turbulent viscosity coefficient (TURB) was varied to 
determine its effect on the fluid motion (Section 2. 1). Figure 17 shows the basic fluid 
motion at approximately 1. 1 sec for three different values of TURB. It is clear that 
the lower the value of TURB the faster is the velocity of the fluid and hence the sooner 
the geyser begins to form. It is emphasized that the turbulent viscosity is calculated in 
a cell containing fluid when at least two of its adjacent neighboring cells also contain 
fluid. This criteria is needed so that 3u/dz and ?>v/9 r can both be computed. Unfort- 
unately due to the coarse grid mesh this condition is not always satisfied when in reality 
it should be. This is true when the fluid is sliding down the wall and is moving across 
the bottom boundary as a layer of fluid with a thickness less than the grid dimension. 
While this condition exists for a relatively short duration during Case 2, this problem 
will become more severe during the remaining cases. 

Since the lack of a surface pressure interpolation scheme has been shown to result in 
surface instability (Figure 15), it was concluded that it is necessary to 
maintain the surface pressure interpolation scheme during geyser formation. Also, 
to minimize surface breakup in the geyser itself it is necessary to add surface marker 
particles over the entire surface to include the geyser. Therefore, a subroutine was 
written to internally add the necessary surface marker particles at selected times. 
Figure 18 shows a case where the subroutine was used to add particles. 

Figures 19A through 19L give computer generated SC4020 plots of the Marker particles 
for Case 1. The surface pressure interpolation scheme was discontinued at 1.566 
seconds. This is the reason that plots 19G through 19L have no solid line outlining 
the surface as is the case in plots 19A through 19F. Figures 20A through 20C give 
computer generated SC4020 plots of the velocity vectors at selected times. 

Figure 21 gives photographs of the drop tower test results simulated by Case 2. 

Initially the interface shape corresponds to a Bond number of 15. The test model was 
a scaled model of the current Centaur tank. Except for the deviation in the initial 
interface shape and the fact that the model had minor differences in tank end shapes, 
the remaining characteristics of settling acceleration, fluid properties, and percent 
liquid are the same as outlined for Case 2. All the times associated with each photo- 
graph in Figure 21 are rounded off to the nearest tenth of a second. 
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A. TURB = 0 B. TURB = .1 C. TURB = 1.0 

Time = 1.025 

Figure 16. Test Model With 30% Ullage Figure 17. Test Model at Approximately 1.1 

With Particle Addition by Subroutine Second for Various Values of TURB 

ADPART 



A. Without Added Surface 
Marker Particles 


B. With Added Surface 
' Marker Particles 


Figure 18. Test Modd. at Approximately 1.375 Seconds With 
and Without Added Surface Marker Particles 
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Figure 19. SC4020 Marker Particle Plots for Case 2 
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Figure 20. SC4020 Velocity Vectors Plots for Case 2 
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Figure 21. Continued 
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The analytically calculated motion was compared with actual drop tower test data 
supplied by NASA/LeRC. The following indicates the close degree of correlation. 

• In both cases the liquid initially moves down the wall of the container with no 
Taylor instability appearing in the center of the tank. The liquid hits the bottom 
of the container at approximately 0. 7 sec (compare Figure 19C with 21B). 

• In both cases the liquid then moved along the bottom surface toward the 
center of the tank. The liquid does not separate from the bottom boundary. 

• When the liquid reached the center of the tank a geyser formed as indi- 
cated in plot 19D. In both cases geyser formation began in less than one 
second (in ... 99 second for the analytical case and approximately .95 
seconds in the test results). 

• The geyser continued to rise and impinged on the fluid remaining in the 
upper portion of the tank. In both cases this impingement occurred in 
approximately 1.25 seconds (see Figure 19E). * In both cases the width 
halfway up the geyser is approximately 3. 0 cm when it impinges the fluid 
in the upper region (see Figure 19E). 

"• In both cases the distance from the top of the tank to the free surface of 
the liquid in the top of the tank is approximately 10 cm when the geyser 
impinges the fluid in the upper region (see Figure 19E). 

The correlation in the detailed fluid motion following geyser impingement . 
on the fluid is not so accurate as above. In the test results the geyser 
actually pushes through the fluid in the upper region and hits the top of 
the tank at approximately 1.45 seconds. As can be seen in plots 19G 
through 19L the geyser does not push through the fluid to reach the top of 
the tank. 

While the detailed motion after 1.25 second cannot be resolved by the program for 
the present grid mesh size it is felt that the lumped motion of the fluid is still accurate 
enough to give a good estimate of the "collected volume" for the duration of this prob- 
lem i Also some other trends are noted. 

• Beginning at about 2 seconds the geyser begins to collapse. This can be 
seen in velocity plot 20C. This collapsing of the geyser causes the liquid 
in the geyser to move as indicated by the arrows in Figure 19J. This 
same type of the geyser breakup is seen in the test results. 

* As can be seen in Figure 19 E there is a space between the geyser and the fluid 
in the upper region of the tank. However, as far as the program is concerned 
impingement has occurred since some of the geyser marker particles are in 
the same cell as the particles defining the surface of the fluid in the upper 
region. 
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• After approximately 1 . 7 seconds (Figure 19H) it becomes increasingly 
difficult to determine the true "collected volume. " By the end of the 
problem at approximately 2.4 seconds there are a number of regions 
which appear to lack Marker particles (see Figure 19L). Does this 
actually imply entrained vapor as can be seen in the test results or is 
it a result of the course grid mesh causing inaccuracies in the fluid 
motion? Recall that once one marker particle gets into a cell the cell 
may be considered full of liquid for purpose of simulation. Also, the 
coarseness of the grid mesh causes the velocity vectors between two 
adjacent cells to deviate by over 90° in some cases. This is particularly 
true when the fluid coming down the tank wall hits the bottom boundary or 
the fluid on the bottom and is forced to change directions (see Figure 20C). 
The result is a dip in the fluid as indicated in Figures 19F and 19K. 

These dips should be neglected when determining the true "collected 
liquid volume. " The lines to the right of some of the SC4020 Marker 
particle plot indicate the height of the liquid in the "collected volume. 

A basic criteria for estimating this height is given in Appendix A. 

Up to approximately 1.416 second (Figure 19F) there is only one line indicating the 
"collected liquid" height. After that there are as many as three lines (if there are 
three the top implies maximum collected volume height while the lower implies the 
minimum collected volume height and the middle one implies the mean or estimated 
actual collected volume height). After seeing the LeRC test results, it is estimated 
that the true "collected liquid" height is between the minimum and the mean height 
while the region between this true height and the maximum height is a region which is 
filled with liquid and entrained vapor. Based on this assumption Figure 22 compares 
the analytical and experimental collected liquid height. It is noticed that the largest 
deviations occur during the initial part of the problem while after 2. 0 seconds the 
analytical and test data agree almost exactly. Due to the ambiguity in any exact 
value for the true height (experimental as well as analytical) it is felt that the program 
does predict the "collected liquid" in this typical settling problem. 

Figure 23 is the collected liquid volume corresponding to the collected liquid height 
in Case 22. 

CASE 3. Case 3 had the following additional characteristics 
20% liquid 

One cell truncated spherical segment bottom dome 
Settling acceleration of 70 cm/sec 2 
Turbulent viscosity coefficient = 0.05 
Initial Marker particle density 7X7 per cell 

Figures 24A through 24F give computer generated SC4020 plots of the Marker particles 
The surface pressure interpolation scheme was discontinued at 0.90 sec (Figure 24E). 
It was found necessary to discontinue this interpolation scheme when the geyser 
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Figure 22. 


Figure 23. 



Comparison of the Analytical and Experimental "Collected 
Liquid"Height for Case 2 


TOTAL LIQUID VOLUME 2401 cm 1 



The Analytically Computed "Collected Liquid Volume" 
Versus Settling Time for Case 2 



impinges on fluid remaining at the top of the tank. Figures 25A and 25G give computer 
generated SC4020 plots of the velocity vectors at selected times. Initially the fluid was 
at rest (all velocities = 0). The scale of the velocity vectors was changed after the 
initial geyser formation because the large velocities in the geyser would go off the plot 
at the initial scale. This is the reason why the velocity vectors of the fluid going down 
the wall are longer in Figure 25A than in 25B. 

The results of this run correlated with actual drop tower test data and existing theory to 
the point of geyser formation. Figure 26 indicates how the leading edge velocity corre- 
lates to existing data (Ref. 2). The reason for the leading edge exhibiting a downward 
acceleration 13% less than freefall is basically because of the noslip condition which 
exists in real fluid flow problems (Ref. 7). However, during all the cases run during 
this contract the freeslip boundary condition exists on all boundaries. It is therefore 
concluded that the fact that the analytical leading edge velocity departs from the settling 
velocity and results in a close agreement with actual leading edge velocity is fortunate. 
The freeslip condition was used instead of the noslip condition since the grid mesh was 
too coarse to resolve the detailed motion of an actual boundary layer. However, the 
coarse grid mesh caused inaccuracies in the numerical computation which resulted in 
the deceleration of the leading edge. It is expected that with more accuracy, resulting 
from a finer grid mesh along the boundary, the leading edge acceleration would approach 
the settling acceleration if the nodip boundary condition continued to be imposed. The 
result would be that the leading edge velocity would be too high. If this occurs in future 
runs it is suggested that the HOPI code be modified so that a partial slip condition can 
be input for cases where the grid mesh is not fine enough for the noslip condition. 

It was noted after running Case 2 that the leading edge velocity was slower than theoreti- 
cally predicted. An examination of this code indicated that some velocities in empty 
cells were being prematurely deleted. This error was corrected before running Case 3. 
It is emphasized that the error was minor for Case 2 and that the basic conclusions 
developed from Case 2 still apply for Case 2. 

The geyser forms initially at . 840 sec. However, the geyser velocity was faster than 
in the drop tower experiments and once the geyser did hit the top of the tank the fluid 
motion could not be accurately resolved. The reason the fluid motion in the 20% liquid 
case cannot be accurately resolved can be seen by comparing Figure 24E with Figure 10. 
At this point the fluid is within thin sheets along the top bulkhead, the side wall, the 
bottom bulkhead and in the geyser. In fact the sheets are thinner than the cell width 
thereby causing almost all the fluid to be within surface cells. This means that the grid 
mesh is not fine enough to resolve the motion. To obtain the necessary resolution to 
accurately solve this problem, with the constant grid mesh restriction that now exists 
in the code, would require at least reducing by half the grid size in both the radial and 
axial direction. This would push the core size of the computer to the limit and would 
more than double the computer time needed to sclve the problem. 

It is therefore concluded that the code in its present form should not be used to solve 
this type of problem. It is recommended that a variable grid mesh option be added to 
the program so that smaller grids can be used where needed (such as around corners). 
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Figure 24. SC4020 Marker Particle Plots for Case 3 



Figure 25. SC4020 Plots of Velocity Vectors for Case 3 
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Figure 26. Correlation of SMAC Leading Edge Velocity With Theory for 
Case 2 (g z = 70 cm/sec 2) 
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Conversely, large grids would be available in the same problem where detailed 
resolution of the motion is not so critical. 

3.2 FULL SCALE TEST RESULTS 

Two analytical runs were made to simulate full-scale Centaur fuel tank conditions. 

Both problems had the following identical characteristics. 

Height to diameter ratio of 2 
Cylindrical tank with a radius of 150 cm 

Fluid (LH 2 ) had a kinematic viscosity v= 1.92 x 10 -3 cm 2 /sec 
Initial interface shape corresponding to an initial Bond number, Bo of 10 
Hemispherical forward dome 

One cell truncated spherical segment bottom dome 
Initially the fluid was at rest (all velocities = 0) 

A grid mesh of DR = DZ = 21.425714 cm similar to that shown in Figure 10 
All boundaries were input with the free slip boundary condition (Section 2.2. 5) 

In all plots presented in this section just half of a tank is shown due to symmetry. The 
left side of each plot corresponds to the centerline of the tank. 

CASE 4. The first full scale run (called Case 4) had the following additional 
characteristics. 

• 70% liquid 

• Settling acceleration of . 27 cm/sec 

• Initial Marker particle density 6x6 per cell (4285 particles) 

Figures 27A through 270 give computer generated SC4020 plots of the Marker particles 
The surface pressure interpolation scheme was discontinued at 50 sec. It has been 
found necessary to discontinue this interpolation scheme when the geyser impinges on 
fluid remaining at the top of the tank. Figures 28A through 28F give computer 
generated SC4020 plots of the velocity vectors at selected times. The scale of the 
velocity vectors was changed after the initial geyser formation because the large 
velocities in the geyser would go off the plot at the initial scale. This explains the 
longer velocity vectors of fluid going down the wall in Figure 28A compared with 28B. 

CASE 5 . The second full scale run (called Case 5) had the following additional 
c harac ter istic s . 

• 85% liquid 

• Settling acceleration of . 135 cm/sec 

• Initial Marker particle density of 5 x 5 per cell (3652 particles) 

Figures 29 A through 29N give computer generated SC4020 plots of the Marker 
particles.' The surface pressure interpolation scheme was discontinued at 55 
seconds. Figures 30A through 30F give computer generated SC4020 plots of the 
velocity vectors at selected times . 
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Figure 29. SC4020 Marker Particle Plots for Case 5 
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3.2.1 DISCUSSION OF RESULTS OF CASE 4. 

• Initially the liquid moves down the wall of the container with a leading edge 
velocity that correlates with existing data (Reference 2). The fluid hits the 
truncated spherical bottom of the tank at 30.5 seconds with a leading edge 
velocity of 7.2 cm/sec. 

• The liquid then travels along the bottom boundary as a thin film. The liquid 
is accelerated as it moves toward the center of the tank so that the velocity 
increases inversely withjhe tank radius. It is noted that turbulent viscosity 
does not yet have any affect at the lower right hand corner. This is because 
the turbulent viscosity (Section 2. 1) is not considered in a corner which is 
surrounded only by surface cells. It is felt that the grid mesh should be finer 
at the lower right hand corner (and all corners) so that the turbulent viscosity 
can take affect. It is suspected that in this analytical run the fluid velocity 
across the bottom boundary is faster than would actually occur if turbulent 
viscosity were presented in the corner. 

• The geyser forms at 45. 3 seconds. Again turbulent effects are lacking because 
of the coarse grid mesh. 

• The geyser forces its way through the fluid in the upper section of the tank as 
indicated in Figures 27C and 27D. The liquid continues to move down the walls 
of the tank and across the bottom boundary. While there is a lack of Marker 
particles next to the bottom tank wall in Figure 27D this region is considered 
to be full of fluid. This is noted in the velocity vectors in Figure 28C. It is 
noted that velocity vectors appear only in cells flagged either FUL or SUR. 

• With increased time all the liquid initially in contact with the upper tank 
boundary is forced down the wall by the up coming geyser (Figures 27E and 
27F). The amount of liquid collected is steadily increasing as indicated by 
Figure 31. 

• The liquid becomes extremely turbulent in the upper region of the tank making 
it unclear if a region should actually be containing entrained vapor or not. 

This is the reason the collected liquid volume is not given after 100 seconds 
in Figure 31. By about 200 seconds (Figure 27M and 27N) the upper center 
region of the tank is clear. However, by 225 seconds (Figure 270) the entire 
tank is considered by the program as full of fluid. This includes even the 
top of the tank. This is because of the coarse grid mesh which causes the 
entrained vapor to be lost in the liquid. This point is discussed in Section 
3.2.3. 
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COLLECTED LIQUID, cm J x 10" 


TOTAL LIQUID VOLUME = 2.47 X 10 7 cm 3 



SETTLING TIME, sec 

Figure 31. Collected Liquid Volume for Case 4 
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It is suggested that after approximately 100 seconds the loss in the calculation 
mesh of entrained bubbles begins to distort the results, since regions that 
should contain vapor are considered to be full of fluid. This distortion in the 
results becomes increasingly worse as more and more of the vapor gets lost 
in the fluid. By 200 seconds the tank is almost completely filled with fluid. At 
this time only the top most point is considered empty. However, there is still 
a sloshing motion of the fluid in the tank. In time (by about 225 seconds), this 
sloshing motion causes the entire tank region to be considered as full of fluid. 

It is noted that after 90 seconds the magnitude of the velocities of the fluid are 
continually decreasing with time. This can be seen by comparing Figures 
28D, E, and F. The magnitude of the velocity vector is directly proportional 
to the length of the velocity vector lines in each plot. At 200 seconds, Figure 
28E, the maximum velocity component in either direction (r or z) is less 
than 4 cm/sec . 

In conclusion, after approximately 150 seconds the results for Case 3 have become 
so distorted by the lack of true empty cells that the results have only marginal 
significance. 

3.2.2 DISCUSSION OF RESULTS FOR CASE 5 . 

v ' • . 

• Initially the motion is the same as for Case 4. The liquid moves down the tank 
wall then accelerates across the bottom boundary. As the liquid converges at 
the center of the tank a geyser forms which impinges on the fluid in the upper 
region of the tank and then forces its way through the fluid. 

• The collected volume (Figure 32) continues to increase as the gas in the tank 
appears to move up in a simple bubble (actually two bubbles since the SC4020 
plots represent only half of a tank). It appears that as the bubble approaches 
the top of tank (until approximately 130 sec) there is not a significant amount of 
break up in the bubble. This implies no entrained vapor in the lower region of 
the tank. 

• By 140 seconds, while it is still felt that the bubble is reasonably intact, the 
shape of the vapor region has been distorted so that actual resolution of the 
gaseous region is not possible because of the coarse grid mesh. This point is 
discussed in the following Section 3.2.3 and in Section 3.2.1. 

3*2.3 GENERAL DISCUSSION OF RESULTS . While r unni ng Case 3 it was concluded 
that the grid mesh used was not fine enough to resolve the fluid motion. This was 
because the cell size was large relative to the liquid available. This resulted in 
nearly all the liquid appearing in cells flagged as SUR cells. While r unni ng Cases 4 
and 5 the same type of problem developed; only in these cases the cell size was large 
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COLLECTED VOLUME, cm a x 10“ 


TOTAL LIQUID VOLUME FOR CASE 4 = 3.06 x 10 7 cm 3 



0 20 40 60 80 100 120 140 

SETTLING TIME, sec 


Figure 32. Collected Liquid Volume for Case 5 
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relative to the gas available. After a given length of time the condition developed 
where there was no empty (G = 4) cells. Figure 33 illustrates how a region containing 
a bubble which should be indicated as partly empty is indicated as a completely full 
section since no single cell is completely empty. 

This is a serious problem since once the entire grid mesh contains no more G = 4 
cells it is impossible to create G = 4 cells. This is because in SMAC a FUL cell 
must be a SUR cell before being empty (flagged G=4) and for a cell to be flagged 
SUR it must have an empty neighbor. In an attempt to remedy this problem, the 
condition was added that a FUL cell could be transformed directly into an EMP cell 
if the full cell lost all its marker particles. While this appears to give more 
meaningful results the abrupt change from FUL to EMP causes significant changes 
m the pressures which greatly increases the number of cycles required for con- 
vergence. Also, the fact remains that a region is at least temporarily flagged as 
full when it is at least partially empty. 

It is concluded that not only should the cell size be small compared to the fluid 
available but should also be small compared with the gas available. This was not 



50 


true for Case 4 which initially had only 13 cells flagged as empty (G = 4) and 7 of 
these G = 4 cells were along the arbitrary boundary which is not even a whole cell. 

Since the problems encountered in Case 2 and Cases 3 and 4 result from the grid 
mesh being too coarse it is concluded that the complex problem of liquid reorientation 
in a tank requires a finer mesh than what has been used during this contract. This 
problem will be even more severe for simulations with internal baffles which are 
characterized by extreme turbulence. It is therefore suggested that HOPI be modified 
to reduce the computer core required to run a problem and that every iteration loop 
be optimized. The core required for problems 1 through 5 was approximately 131, 000 
for the initial setup (with ARBND) and 123,000 for successive runs (without ARBND 
and VOLUME). It is felt that by using "overlay" the core could be reduced to below 
100,000. This would not only reduce computer cost but would also make available 
the additional core required when a finer grid mesh is needed. 



4.0 CONCLUSIONS AND REC OMME NDATIONS 


During this contract a computer program called HOPI was developed. This program 
uses the Simplified Marker and Cell, SMAC, numerical technique. HOPI has shown 
that it can analytically determine the fluid motion in a cryogen storage tank under a 
continuous settling load. The equations programmed in HOPI are applicable to 
Newtonian incompressible fluid. HOPI can be used with curved boundaries. While 
the grid dimension can differ in the radial, DR, and axial, DZ, direction; HOPI can 
presently handle only a constant-size grid mesh. Experience has indicated that this 
is a serious limitation since a finer grid mesh is usually needed at corners and along 
the boundaries than is needed in the middle region of the tank. It would be extremely 
inefficient to use the smallest needed cell throughout the grid mesh. Currently HOPI 
and the related subroutines are loaded into the computer core at the same time. This 
causes an unnecessarily high demand for core. In fact, the grid size used for a 
typical problem loads the core to near capacity. It is estimated that if the radial 
and grid dimension were both reduced by half, thereby multiplying by a factor of four 
the grid size, the core of the computer would be exceeded. However, experience 
indicated that a reduction of the grid dimensions by at least half is needed to accurately 
resolve the fluid motion under a number of conditions. This is especially true if 
there is less than 25% liquid or gas or if the liquid becomes extremely turbulent. 

HOPI has a surface pressure interpolation scheme which is required to avoid 
unrealistic surface breakup. This scheme sets the pressure in the middle of a 
surface cell so that the pressure at the actual surface of the liquid is correct. 

HOPI uses the equation 

Vb =TroBx ^ max (iHi- i»4-i) 

to simulate a turbulent viscosity, where 



and 

v is the radial component of velocity 
u is the axial component of velocity 
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TURB is an empirical input quality . It has been found that the use of a proper 
value of TURB does aid in correlating the analytical results with the test results. 
Experience has also indicated that Hie determination of the "collected liquid 
volume" is best determined from an examination of the SC4020 Marker particle 
plots. Lastly, the subroutine written to determine "collected liquid volume" has 
not proven very useful, especially where turbulence causes vapor entrainment in 
the liquid. 

For the future work, the following changes in HOP I are recommended. 

1. The subscripting should be changed from a double to a single throughout the 
code, This will reduce the computer running time for a given problem. 

2. Use overlay to reduce the core needed for a given size problem. 

3. Delete subroutine VOLUME. Because of its limited usefulness, it appears that 
it can be deleted. Furthermore this will reduce the core required for a given 
problem. 

4. For a half tank simulation, set-up the case with finer grid than the 288 cells used 
to define the tank during this contract. 

5. Introduce a variable grid mesh. This will allow more efficient use of the 
computer by using smaller grid sizes where needed but allowing larger size 
grids where possible. 

6. Develop a surface pressure interpolation scheme which can be used even after 
the geyser impinges on the liquid in the upper region of the tank. 
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APPENDIX A 


SELECTION CRITERIA FOR THE MINIMUM, MEAN AND 
MAXIMUM COLLECTED LIQUID HEIGHT 

The definition of the "collected liquid " height is the height below which liquid and 
no vapor exist. As can be seen in the Marker particle plots this height is often 
difficult to define. This is because a cell is considered full once it contains a single 
marker particle. Therefore, the fact that the region looks nearly free of particles 
or at least the particles are rarified when compared to other regions does not in 
itself indicate that vapor exists in this region. Secondly, the coarse grid mesh, 
while needed to avoid excessive time steps during computation, results in inaccur- 
acies in the fluid motion. When the velocity vectors of adjacent cells deviate by 
more than 90° (see Figure 20C) the result is an unrealistic dip in the fluid as indi- 
cated in Figure 19H. 

A basic criteria used to estimate the collected liquid height is presented in the following 
paragraphs. It is used only as a guide and is not an absolute criteria. 

Up to 1.416 seconds (Figure 19F) there is only one height. Initially this height was 
computed using the assumption that the height equalled the actual minimum level of 
the particles. Then as the dip in the surface became more pronounced as in Figure 
19F, the dip portion was considered as actually containing particles. 

Starting with Figure 19G more than one height is indicated because of ambiguity in 
what is the true "collected liquid" height. Starting with Figure 19H three heights are 
usually indicated. These heights were estimated as follows. 

Minimum Height. Assumes that the dominant lowest dip in the surface is actually 
filled with particles to the point where there is an abrupt change in direction in the 
particles outlining the surface of the geyser. In Figure 19H this point is noted by 
( 1 ). 

Mean Height. Assumes that particles actually fill up the void to the point where the 
particles nearest the bottom yet in the fluid coming down the wall start to move to- 
ward the center of the tank instead of directly toward the bottom. In Figure 19H this 
point is noted by (2). Or, as is the case in Figure 19L, the mean height is estimated 
as existing at the highest point of the fluid initially from the side fluid and which was 
not part of the geyser fluid. / 
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Maximum Height . There are two methods for estimating the maximum height. One 
assumes that the "collected liquid" is to the point where the particles in the fluid 
coming down the wall first start to move toward the center of the tank instead of 
directly down. In Figure 19H this point is noted by (3). The other assumes that 
the height exists where the geyser has collapsed and flowed toward the liquid on 
the side of the tank and is in contact with the fluid coming down the side of the tank 
wall. Contact exists when particles of the geyser side and particles in the wall 
side occupy part of the same cell. This criteria is illustrated in Figure 19 J. 
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NOMENCLATURE 


area 

Bond number, gR 2 /yS 

velocity divergence = (l/r a ) OrV^r) + (Sv/az) 
grid dimension in the radial direction 
grid dimension in the axial direction 
cell flag 

radial acceleration 

axial acceleration 

radial spatial coordinate index 

axial spatial coordinate index 

boundary slope 

time index number of cycles 

unit normal defining a boundary segment 

pressure 

pressure of cell ij 
radial coordinate 
tank radius 

relaxation parameter 
Poisson source term 
empirical factor for turbulence 
time 

radial component of velocity 

radial storage variable, radial component of the tilde velocity 
axial component of velocity 

axial storage variable, axial component of the tilde velocity 
velocity 





liquid velocity at midpoint of boundary segment computed with 
the MAC area weighing interpolation scheme . 

position of midpoint of a boundary segment 
position of particle 


% axial coordinate 

a geometric parameter, a = 1. 0 in cylindrical coordinates and 

equals 0.0 in plane (cartesian) coordinates 

/? specific surface tension 

6 r incremental step in the r direction 

6 z incremental step in the z direction 

At time step 

Ar radial mesh width 

Az axial mesh width 

€ boundary sensing parameter 

X minimum mesh dimension, minimum of DR or DZ 

v kinematic viscosity 

p density 

0 true pressure normalized to unit density 

a At/p 

i]j arbitrary pressure normalized to unit density (pseudopressure) 

oo vorticity 

v gradient operator 

Superscripts 

k iteration index 

n counts time cycles 

Subscripts 


i position in the finite-difference mesh 

j position in the finite-difference mesh 

p located at midpoint of boundary segment 
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